R Scott Collings
Math 321 Sec 002
27 Nov 2018
from matplotlib import pyplot as plt
from scipy.io import wavfile
from scipy import fftpack
import numpy as np
import IPython
from copy import deepcopy
from scipy import signal
from imageio import imread
plt.rcParams["figure.dpi"] = 300 # Fix plot quality.
plt.rcParams["figure.figsize"] = (12,3) # Change plot size / aspect (you may adjust this).
class SoundWave(object):
"""A class for working with digital audio signals."""
# Problem 1.1
def __init__(self, rate, samples):
"""Set the SoundWave class attributes.
Parameters:
rate (int): The sample rate of the sound.
samples ((n,) ndarray): NumPy array of samples.
"""
self.rate = rate
self.samples = samples
# Problems 1.1 and 1.7
def plot(self, transform=False):
"""Plot the graph of the sound wave (time versus amplitude)."""
x = np.linspace(0,len(self.samples)/self.rate,len(self.samples))
if (transform):
plt.subplot(121)
plt.ylim(-32768,32767)
plt.xlabel('sample index')
plt.ylabel('sample')
plt.gca().set_title('Audio File')
plt.plot(x,self.samples)
if(transform):
dft = np.abs(fftpack.fft(self.samples))
plt.subplot(122)
plt.xlim(0,self.rate / 2)
plt.plot(x * self.rate * self.rate / len(self.samples),dft)
plt.gca().set_title('Audio File DFT')
plt.xlabel('frequency')
plt.ylabel(r'$c_k$')
plt.show()
# Problem 1.2
def export(self, filename, force=False):
"""Generate a wav file from the sample rate and samples.
If the array of samples is not of type np.int16, scale it before exporting.
Parameters:
filename (str): The name of the wav file to export the sound to.
"""
samples = deepcopy(self.samples)
if(type(self.samples[0]) != np.int16 or force):
samples = np.real(samples * 32768.0 / np.max(np.abs(samples))).astype(np.int16)
wavfile.write(filename,self.rate,samples)
# Problem 1.4
def __add__(self, other):
"""Combine the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to add
to the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the combined samples.
Raises:
ValueError: if the two sample arrays are not the same length.
"""
if(len(self.samples) != len(other.samples) or self.rate != other.rate):
raise ValueError("Sound files must be same length and sample rate")
return SoundWave(self.rate,self.samples + other.samples)
# Problem 1.4
def __rshift__(self, other):
"""Concatentate the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to concatenate
to the samples contained in this object.
Raises:
ValueError: if the two sample rates are not equal.
"""
if(self.rate != other.rate):
raise ValueError('Sound files must have same sample rate')
return SoundWave(self.rate,np.append(self.samples,other.samples))
# Problem 2.1
def __mul__(self, other):
"""Convolve the samples from two SoundWave objects using circular convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
if (self.rate != other.rate):
raise ValueError('Sample rates must match')
diff = len(self.samples) - len(other.samples)
samples = []
f_hat = []
g_hat = []
if (diff > 0):
samples = np.append(deepcopy(other.samples), np.zeros(diff))
f_hat = fftpack.fft(self.samples)
g_hat = fftpack.fft(samples)
else:
samples = np.append(deepcopy(self.samples), np.zeros(abs(diff)))
f_hat = fftpack.fft(samples)
g_hat = fftpack.fft(other.samples)
return SoundWave(self.rate,fftpack.ifft(np.multiply(f_hat,g_hat)))
# Problem 2.2
def __pow__(self, other):
"""Convolve the samples from two SoundWave objects using linear convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
if (self.rate != other.rate):
raise ValueError('Sample rates must match')
n = len(self.samples)
m = len(other.samples)
length = 1
while length < n + m - 1:
length *= 2
samples1 = SoundWave(self.rate,np.append(self.samples,np.zeros(length - n)))
samples2 = SoundWave(other.rate,np.append(other.samples,np.zeros(length - m)))
SW = samples1 * samples2
return SoundWave(SW.rate,SW.samples[:n + m - 1])
# Problem 2.4
def clean(self, low_freq, high_freq):
"""Remove a range of frequencies from the samples using the DFT.
Parameters:
low_freq (float): Lower bound of the frequency range to zero out.
high_freq (float): Higher boound of the frequency range to zero out.
"""
low_freq = low_freq * len(self.samples) // self.rate
high_freq = high_freq * len(self.samples) // self.rate
dft = fftpack.fft(self.samples)
self.samples = fftpack.ifft(
np.append(
np.append(dft[:low_freq],np.zeros_like(dft[low_freq:high_freq])),
np.append(
np.append(dft[high_freq:-high_freq],np.zeros_like(dft[-high_freq:-low_freq])),
dft[-low_freq:])
)
)
SoundWave.__init__().SoundWave.plot().scipy.io.wavefile.read() and the SoundWave class to plot tada.wav.rate,samples = wavfile.read('tada.wav')
SW = SoundWave(rate,samples)
SW.plot()
SoundWave.export().export() method to create two new files containing the same sound as tada.wav: one without scaling, and one with scaling (use force=True).IPython.display.Audio() to embed the original and new versions of tada.wav in the notebook.rate,samples = wavfile.read('tada.wav')
SW = SoundWave(rate,samples)
SW.export('tadaNotScaled.wav')
SW.export('tadaScaled.wav',force=True)
IPython.display.Audio('tada.wav')
IPython.display.Audio('tadaNotScaled.wav')
IPython.display.Audio('tadaScaled.wav')
generate_note().generate_note() to create an A tone that lasts for two seconds. Embed it in the notebook.def generate_note(frequency, duration):
"""Generate an instance of the SoundWave class corresponding to
the desired soundwave. Uses sample rate of 44100 Hz.
Parameters:
frequency (float): The frequency of the desired sound.
duration (float): The length of the desired sound in seconds.
Returns:
sound (SoundWave): An instance of the SoundWave class.
"""
x = np.linspace(0,duration,44100 * duration)
f = np.sin(2 * np.pi * frequency * x)
return SoundWave(44100,f)
SW = generate_note(440,2)
SW.export('a4note.wav')
IPython.display.Audio('a4note.wav')
SoundWave.__add__().SoundWave.__rshift__().SWA = generate_note(440,3)
SWC = generate_note(523.25,3)
SWE = generate_note(659.25,3)
SW = SWA + SWC + SWE
SW.export('aceChord.wav')
IPython.display.Audio('aceChord.wav')
SWA = generate_note(440,1)
SWC = generate_note(523.25,1)
SWE = generate_note(659.25,1)
SW = SWA >> SWC >> SWE
SW.export('aceArpeggio.wav')
IPython.display.Audio('aceArpeggio.wav')
simple_dft() with the formula for $c_k$ given below.np.allclose() to check that simple_dft() and scipy.fftpack.fft() give the same result (after scaling).def simple_dft(samples):
"""Compute the DFT of an array of samples.
Parameters:
samples ((n,) ndarray): an array of samples.
Returns:
((n,) ndarray): The DFT of the given array.
"""
n = len(samples)
x = range(n)
W = (1/n) * np.array([[np.exp(-2 * np.pi * 1j * l * k / n) for l in x] for k in x])
return W @ samples
%time mine = simple_dft(np.array([1,6,2,4]))
print(mine)
theirs = fftpack.fft(np.array([1,6,2,4]))
print(theirs)
np.allclose(mine, theirs / len(theirs))
simple_fft().simple_dft(), simple_fft(), and scipy.fftpack.fft().
Print the runtimes of each computation.np.allclose() to check that simple_fft() and scipy.fftpack.fft() give the same result (after scaling).def simple_fft(samples, threshold=1):
"""Compute the DFT using the FFT algorithm.
Parameters:
samples ((n,) ndarray): an array of samples.
threshold (int): when a subarray of samples has fewer
elements than this integer, use simple_dft() to
compute the DFT of that subarray.
Returns:
((n,) ndarray): The DFT of the given array.
"""
def split(g):
n = len(g)
if (n <= threshold):
return simple_dft(g)
else:
even = split(g[::2])
odd = split(g[1::2])
z = np.array([np.exp(-2 * np.pi * 1j * y / n) for y in range(0,n)])
m = n // 2
return np.append((even + np.multiply(z[:m], odd)), (even + np.multiply(z[m:], odd)))
return split(samples) / len(samples)
x = np.random.rand(8192)
%time bad = simple_dft(x)
%time mine = simple_fft(x)
%time theirs = fftpack.fft(x)
np.allclose(mine * len(mine),theirs)
SoundWave.plot() so that it accepts a boolean. When the boolean is True, take the DFT of the stored samples and plot (in a new subplot) the frequencies present on the $x$-axis and the magnituds of those frequences on the $y$-axis. Only the display the first half of the plot, and adjust the $x$-axis so that it correctly shows the frequencies in Hertz.rate,samples = wavfile.read('a4note.wav')
SW = SoundWave(rate,samples)
SW.plot(transform=True)
rate,samples = wavfile.read('aceChord.wav')
SW = SoundWave(rate,samples)
SW.plot(transform=True)
Use the DFT to determine the individual notes that are present in mystery_chord.wav.
rate,samples = wavfile.read('mystery_chord.wav')
SW = SoundWave(rate,samples)
SW.plot(transform=True)
x = fftpack.fft(SW.samples)
notes = np.where(np.abs(x[:len(SW.samples) // 2]) > 5e8)[0]
notes * SW.rate / len(SW.samples)
The notes are...A C D G
SoundWave.__mul__() for circular convolution.tada.wav.tada.wav and the white noise. Embed the result in the notebook.rate,samples = wavfile.read('tada.wav')
SWTada = SoundWave(rate,samples)
whiteNoise = SoundWave(rate,np.random.randint(-32767, 32767, rate * 2, dtype=np.int16))
SW = SWTada * whiteNoise
SW.export('whiteNoiseTada.wav')
IPython.display.Audio('whiteNoiseTada.wav')
SWLoop = SW >> SW
SWLoop.export('loopedWhiteNoiseTada.wav')
IPython.display.Audio('loopedWhiteNoiseTada.wav')
SoundWave.__pow__() for linear convolution.CGC.wav and GCG.wav using SoundWave.__pow__() and scipy.signal.fftconvolve().SoundWave.__pow__() and scipy.signal.fftconvolve() sound the same.rate1,samples1 = wavfile.read('CGC.wav')
rate2,samples2 = wavfile.read('GCG.wav')
SWCGC = SoundWave(rate1,samples1)
SWGCG = SoundWave(rate2,samples2)
%time SW1 = SWCGC ** SWGCG
print(samples1)
print(samples2)
%time convolved = signal.fftconvolve(samples1,samples2)
SW2 = SoundWave(rate1,convolved)
SW1.export('powLinConvolve.wav')
SW2.export('scipyLinConvolve.wav')
IPython.display.Audio('CGC.wav')
IPython.display.Audio('GCG.wav')
IPython.display.Audio('powLinConvolve.wav')
IPython.display.Audio('scipyLinConvolve.wav')
Use SoundWave.__pow__() or scipy.signal.fftconvolve() to compute the linear convolution of chopin.wav and balloon.wav.
Embed the two original sounds and their convolution in the notebook.
IPython.display.Audio('chopin.wav')
IPython.display.Audio('balloon.wav')
rate1,samples1 = wavfile.read('chopin.wav')
rate2,samples2 = wavfile.read('balloon.wav')
samples = signal.fftconvolve(samples1,samples2)
SW = SoundWave(rate1,samples)
SW.export('chopinBalloon.wav')
IPython.display.Audio('chopinBalloon.wav')
SoundWave.clean().noisy1.wav by filtering out frequencies from $1250$-$2600$ Hz. Embed the original and the cleaned versions in the notebook.noisy2.wav. Embed the original and the cleaned versions in the notebook.IPython.display.Audio('noisy1.wav')
rate,samples = wavfile.read('noisy1.wav')
SW = SoundWave(rate,samples)
SW.plot(transform=True)
SW.clean(1250,2600)
SW.plot(transform=True)
SW.export('clean1.wav')
IPython.display.Audio('clean1.wav')
IPython.display.Audio('noisy2.wav')
rate,samples = wavfile.read('noisy2.wav')
SW = SoundWave(rate,samples)
SW.plot(transform=True)
SW.clean(1250,4500)
SW.plot(transform=True)
SW.export('clean2.wav')
IPython.display.Audio('clean2.wav')
vuvuzela.wav by filtering bad frequencies out of the left and right channels individually.rate,samples = wavfile.read('vuvuzela.wav')
samples1 = samples.T[0]
samples2 = samples.T[1]
SW1 = SoundWave(rate,samples1)
SW1.plot(transform=True)
SW2 = SoundWave(rate,samples2)
SW2.plot(transform=True)
SW1.clean(200,500)
SW1.plot(transform=True)
SW2.clean(200,500)
SW2.plot(transform=True)
wavfile.write('noVuvuzela.wav',rate,np.array([SW1.samples,SW2.samples]).astype(np.int16).T)
IPython.display.Audio('noVuvuzela.wav')
license_plate.png so that the year printed on the sticker in the bottom right corner of the plate is legible.#im = imread('noisy_face.png')
#print("I can't find license_plate.png anywhere")
#plt.imshow(im,cmap='gray')
#plt.show()
The year on the sticker is...